subroutine rhs(u,u_rhs)
implicit none
real*8::u(:),u_rhs(:)
integer::i,n
real*8,parameter::c=1.,dx=1.

n=size(u)
do i= 4, size(u)-1
    
   u_rhs(i)=-1./12.*c*(3.*u(i+1)+10.*u(i)-18.*u(i-1)+6.*u(i-2)-u(i-3))/(dx**2)
end do

u_rhs(1)=-1./12.*c*(3.*u(2)+10.*u(1)-18.*u(n)+6.*u(n-1)-u(n-2))/(dx**2)
u_rhs(2)=-1./12.*c*(3.*u(3)+10.*u(2)-18.*u(1)+6.*u(n)-u(n-1))/(dx**2)
u_rhs(3)=-1./12.*c*(3.*u(4)+10.*u(3)-18.*u(2)+6.*u(1)-u(n))/(dx**2)
u_rhs(n)=-1./12.*c*(3.*u(1)+10.*u(n)-18.*u(n-1)+6.*u(n-2)-u(n-3))/(dx**2)
end subroutine

program main
implicit none
integer,parameter::n=100,m=500
real*8:: u(n,m)=0,u_k1(n)=0,u_k2(n)=0,u_k3(n)=0,u_k4(n)=0,u_tmp(n)=0
integer::i,j
real,parameter::dt=0.25

interface
  subroutine rhs(u,u_rhs)
  implicit none
  real*8::u(:),u_rhs(:)
  integer::i 
  end subroutine

end interface

!initial data
do i=1,40

  u(i,1)=1

end do

!runge kutta
do i=1,m-1

 call rhs(u(:,i),u_k1)
  do j=1,n
    u_tmp(j)=u(j,i)+0.5*dt*u_k1(j)
  end do 

  call rhs(u_tmp(:),u_k2)
  do j=1,n
    u_tmp(j)=u(j,i)+0.5*dt*u_k2(j)
  end do

  call rhs(u_tmp(:),u_k3)
  do j=1,n
    u_tmp(j)=u(j,i)+dt*u_k3(j)
  end do

  call rhs(u_tmp(:),u_k4)

  do j=1,n
     u(j,i+1)=u(j,i)+(1./6.)*dt*(u_k1(j)+2*u_k2(j)+2*u_k3(j)+u_k4(j))
  end do

!print*,u_k1(41),u_k2(41),u_k3(41),u_k4(41)
!print*,u_k1(1),u_k2(1),u_k3(1),u_k4(1)

end do

open(file='a.txt',unit=10,action="write")
do i=1,m
 write(10,"(100f9.5)")u(:,i)
end do

end program
